/*=========================================================================
 Copyright (c) 2008, Nokia Corporation and/or its subsidiary(-ies).
 All rights reserved.

 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions are met:
 
     * Redistributions of source code must retain the above copyright
       notice, this list of conditions and the following disclaimer.
     * Redistributions in binary form must reproduce the above copyright
       notice, this list of conditions and the following disclaimer in the
       documentation and/or other materials provided with the distribution.
     * Neither the name of the Nokia Corporation and/or its subsidiary(-ies)nor the
       names of its contributors may be used to endorse or promote products
       derived from this software without specific prior written permission.

 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
 EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
 MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
 THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
 OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR 
 TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=========================================================================================*/

#ifndef EIGENVALUES_H
#define EIGENVALUES_H
#include <e32base.h>
#include "ncvMatrix.h"
#include "ncvFixed.h"
/**
 * \class CEigenvalues
 * \brief Eigenvalue decomposition of a real matrix.
 *
 * Usage:
 * <pre>
 * CNokiaCVMatrix * matrix = ...;
 * CEigenvalues * eigs;
 * TRAPD(err, eigs = CEigenvalues->NewL(matrix));
 * if(err != KErrNone)
 *    {
 *    //eigs not found do something
 *    }
 * else
 *    {
 *    TReal * real = eigs->RealEigenvalues();
 *    delete eigs;
 *    }
 * ...
 *
 * </pre>
 */
class CEigenvalues: public CBase
	{

private:
	CEigenvalues();

	/*!
	 * \brief Construct the eigenvalue decomposition of the given matrix
	 * \note the function can leave (ConstructL)
	 */
	void ConstructL(CNokiaCVMatrix &aMatrix);
	/*!
	 * \brief Generate the eigenvalues into the iD and iE arrays (iD real values,
	 *        iE the corresponding imaginary values)
	 *
	 * Eigenvector generation is switched off and not fully implemented,
	 * because there is no support for complex eigenvectors
	 * \note the function can leave
	 */
	void Decompose();


	/*!
	 * \brief Generate the householder reduction of the real SYMMETRIC matrix given in
	 *        the constructor.
	 *
	 * The decomposition algorithm is derived from the Tred2 algorithm from
	 * Wilkinson J.H & Reinsch C, Handbook
	 * for Automatic computing Vol II, Linear Algebra
	 */
	void Tred2()
	;
	/*!
	 * \brief QL algorithm with implicit shift. Finds the eigenvalues (and vectors)
	 *        of the tridiagonal matrix generated by Tred2().
	 *
	 * The decomposition algorithm is derived from the Tqli algorithm from
	 * Wilkinson J.H & Reinsch C, Handbook
	 * for Automatic computing Vol II, Linear Algebra
	 * \note the function can leave
	 */
	void Tqli();

	/*!
	 * \brief Finds the eigenvalues (and vectors, not complex) of a real matrix with
	 *        QR algorithm. The matrix has to be reduced into upper Hessenberg form
	 *        With the methods Ortans and Orthes

	 * \note KNokiaCVErrToManyIterations thrown (ncvErrors.h), if eigs not found fast enough.
	 *
	 * The decomposition algorithm is derived from the Hqr algorithm from
	 * Wilkinson J.H & Reinsch C, Handbook for Automatic computing Vol II,
	 * Linear Algebra.pages 383 - 388 enhanced with
	 * Matlab's new adhoc shift according to JAMA, the java matrix package.
	 *
	 * \note the function can leave
	 *
	 */
	void Hqr2();

	/*!
	 * \brief Generates an unsymmetric matrix into Hessenberg form with orthogonal transforms
	 *        together with the method Ortrans() for the Hqr2 algorithm to use.
	 *
	 *
	 *  \note KNokiaCVErrToManyIterations thrown (ncvErrors.h), if eigs not found fast enough.
	 *
	 *  Derived form the algol procedure Orthes from Wilkinson J.H & Reinsch C, Handbook
	 *  for Automatic computing Vol II, Linear Algebra, pages 349-350
	 *
	 *
	 */
	void Orthes();

	/*!
	 * \brief Forms the accumulated transformations from the info. left from Orthes.
	 *
	 *  Derived form the algol procedure ortans from Wilkinson J.H & Reinsch C, Handbook
	 *  for Automatic computing Vol II, Linear Algebra, page 382
	 */
	void Ortrans();

	/*!
	 * \brief Complex division of the given aX:s., used by Hqr2()
	 *
	 *  Derived form the algol procedure cdiv from Wilkinson J.H & Reinsch C, Handbook
	 *  for Automatic computing Vol II, Linear Algebra.
	 */
	void Cdiv(TReal aXr,TReal aXi, TReal aYr, TReal aYi, TReal *aResR, TReal *aResI);

	TReal Pythag(TReal a, TReal b);

	/*!
	 * \brief Getter for the eigenvectors, if in use.
	 *
	 * If vectors complex only the other one (complement) is added into the
	 * matrix so that the real part is stored into the column i and
	 * the complex part into the column i + 1
	 */
	CNokiaCVMatrix* Matrix();
public:

	IMPORT_C ~CEigenvalues();
	/*!
	 * \brief Construct the eigenvalue decomposition of the given matrix
	 */
	IMPORT_C static CEigenvalues* NewL(CNokiaCVMatrix &aMatrix);

	/*!
	 * \brief Getter for the real part of the eigenvalues (iD)
	 * \return the real part of the eigenvalues in an array
	 */
	IMPORT_C TReal* RealEigenvalues();
	/*!
	 * \brief Getter for the imaginary part of the eigenvalues (iE)
	 * \return the imaginary part of the eigenvalues in an array
	 */
	IMPORT_C TReal* ImaginaryEigenvalues();


private:
	TInt iSize;
	TReal* iD;
	TReal* iE;
	TBool iSymmetric;
	TReal iEps;
	TReal** iH; //Eigenvectors will be stored here
	TReal iTol;
	TReal** iV; //Storage for accumulated transformations with unsymmetric case
	TBool iVectors; //A flag that indicates whether or not to get eigenvectors, this is
	//set to EFalse, due to the fact that complex eigenvectors do
	//not produce correct results and no complex representation implemented
	}
;

#endif // CEigenvalues_H
